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ABSTRACT 

Solar-type stars exhibit a rich variety of magnetic activity. Seeking to explore the convective origins of 
this activity, we have carried out a series of global 3D magnctohydrodynamic (MHD) simulations with 
the anelastic spherical harmonic (ASH) code. Here we report on the dynamo mechanisms achieved 
as the effects of artificial diffusion are systematically decreased. The simulations are carried out at 
a nominal rotation rate of three times the solar value (SfJ©), but similar dynamics may also apply 
to the Sun. Our previous simulations demonstrated that convective dynamos can build persistent 
toroidal flux structures (magnetic wreaths) in the midst of a turbulent convection zone and that high 
rotation rates promote the cyclic reversal of these wreaths. Here we demonstrate that magnetic cycles 
can also be achieved by reducing the diffusion, thus increasing the Reynolds and magnetic Reynolds 
numbers. In these more turbulent models, diffusive processes no longer play a significant role in the 
key dynamical balances that establish and maintain the differential rotation and magnetic wreaths. 
Magnetic reversals are attributed to an imbalance in the poloidal magnetic induction by convective 
motions that is stabilized at higher diffusion levels. Additionally, the enhanced levels of turbulence 
lead to greater intermittency in the toroidal magnetic wreaths, promoting the generation of buoyant 
magnetic loops that rise from the deep interior to the upper regions of our simulated domain. The 
implications of such turbulence-induced magnetic buoyancy for solar and stellar flux emergence are 
also discussed. 

Keywords: starsrinteriors - Sununterior 



> 1 1. MAGNETIC VARIABILITY IN THE SUN AND IN 

; SUN-LIKE STARS 

Informed and motivated by advances in observations of 
,— ' \ the Sun and solar analogues as well as the extensive the- 
. oretical framework of convective dynamo theory, we have 
undertaken a scries of three-dimensional magnctohydro- 
dynamic (MHD) simulations of convection and dynamo 
fSJ ■ action in solar-type stars. These numerical experiments 
t-H ■ show that a rich variety of temporal variability in the 
IL/ 1 magnetic topology, polarity, and field strength can be 
. £h . achieved without varying the rotation rate, even over 
■ the portion of parameter space accessible to our com- 
^ " putational resources. We focus on models with rotation 
C$ \ rates faster than the current solar rate Q@, which accen- 
tuates the magnetic self-organization processes we are 
interested in exploring and which taps into the abun- 
dant observations of magnetic activity in young, rapidly- 
rotating, solar-like stars. 

Our work builds upon three previous groups of simu- 
lations that show the same stellar configuration, namely 
considering dynamics within the deep convective enve- 
lope of our current su n, but having thes e nominally 
young stars rotate faster. iBrown et al.l (|2008f) began with 
hydrodynamic simulations involving a range of rotation 
rates up to 10f2©, finding that strong differential rota- 
tion is realized, and that the columnar convection at 
low latitudes can exhibit significant modulation in am- 
plitude with longitude, even appearing as nearly isolated 
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active nests. IBrown et al.l (|2010() examined dynamo ac- 
tion achieved in a MHD simulation carried out at 3f2©, 
finding that the convection can build ordered global-scale 
magnetic fields that appear as two wreaths of strong 
toroidal field positioned above and below the equator. 
These striking structures can persist for long intervals 
despite being embedded within a turbulent convective 
layer. Turning to dyn amo action proceedi ng at a faster 
rotation rate of 5f2©, IBrown et al.l (1201 If ) showed that 
sclf-consistcntly generated magnetic wreaths at low lat- 
itudes can undergo reversals in global magnetic polarity 
and even quasi-cycles of magnetic activity. The complex 
steps involved in the magnetic field reversals are accom- 
panied by variations in the differential rotation, including 
bands of relatively fast and slow fluid propagating toward 
the poles. 

As we decrease the dissipation in our simulations, we 
find that reversals of global magnetic polarity and cycles 
of magnetic activity are achieved. Despite more vigor- 
ous small-scale turbulence, these simulations still form 
global-scale magnetic wreaths in the bulk of their con- 
vective layers. Yet, the dynamical balances which main- 
tain differential rotation and generate the mean toroidal 
magnetic field change as resolved turbulent transport as- 
sumes the dissipative role that was previously played by 
unresolved subgrid-scale (SGS) turbulent diffusion. We 
attribute the origin of magnetic cycles to a similar shift 
in the dynamical balance in the mean poloidal induction 
equation that had previously sustained steady wreaths. 
As the SGS dissipation is decreased, it is unable to off- 
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set the zonal component of turbulent electromotive force 
(EMF), which generates opposing poloidal field near the 
equator and thereby brings about the polarity reversal of 
the wreaths. A decrease in SGS dissipation also makes 
the wreaths more localized as coherent wreath segments 
over a limited range of longitudes as opposed to axisym- 
metric bands. This tends to decrease the mean field com- 
ponent while simultaneously increasing the field strength 
in the core of the wreaths. We argue that this has im- 
portant implications for flux emergence, since it is these 
localized wreath cores that are most likely to trigger the 
magnetic buoyancy instability. 

1.1. Inspiration from Observational Advances 

As these simulations demonstrate, a rich variety of 
temporal variability can be achieved in dynamo mod- 
els with only modest changes in their control parame- 
ters. Magnetic activity is a ubiquitous characteristic of 
sun-like stars and many stars exhibit cycles of magnetic 
activity. The best example is perhaps the sun's 22-year 
magnetic activity cycle. The interplay of turbulent con- 
vection, rotation, and stratification in the solar convec- 
tion zone creates a cyclic dynamo which drives variations 
in the interior, on the surface, and throu ghout the sun's 
extended atmosphere l|Pinto et al.ll201 lT ) . 

The Sun is not alone in its magnetic variability. Solar- 
type stars generate magnetism almost without exception, 
particularly at rotation rates greater than that of the cur- 
rent sun. Young, rapidly rotating suns appear to have 
much stronger magnetic fields at their surfaces. Obser- 
vations reveal a clear correlation between rotation and 
magnetic activity, as inferred from proxies such as X-ray 
and chromospheric emission dSaar fe Brandenburgj[l999l; 
iPizzolato et al.ll2003t I Wright et al.ll201lf r however super- 
imposed on this trend is considerable variation in the 
presence and period of magnetic activity cycles. There 
have been a number of attempts to monitor the mag- 
netic activity cycles of other stars using solar-calibrated 
proxies for magneti c acti vity (e.g.. iBaliunas et al.lll995t 
iHempelmann et alJll996l : lOlah et al.ll2009f ). These pro- 
grams require long, sustained periods of consistent ob- 
servations, and are therefore rare. To date, the largest 
such project is the Mount Wilson HK survey, which mea- 
sured chromospheric calcium lines as a proxy for mag- 
netic activity for 111 solar- like stars over a 25-year pe- 
riod ending in 1991. In that study almost half of the 
stars showed cyclic behavior incl uding 21 with regula r 
periods between 7 and 25 years (jBaliunas et alJ[l995[ ). 
The existence of sun-like stars without clear cycles of 
magnetic activity provide inspiration in this work for 
the study of a family of dynamo models that lie very 
closely together in parameter space but exhibit markedly 
different degrees of temporal variability in their large- 
scale magnetic features. Improved observational tech- 
niques including spot - tracking from Kep ler photometry 
(jMeibom et al J 1201 It iLlama et all 120121) and Zeeman- 
Doppler imaging (iPetit et al. M200a IGaulme et al.ll20"lH: 
iMorgenthaler et all 120121 )" may permit measurements of 
the size, frequency, and magnetic flux of starspots and 
the topology and spatial variability of photospheric mag- 
netic fields. These developments are likely to provide new 
challenges to our existing theories of dynamo action and 
flux emergence in sun-like stars. 

Whatever the future of solar and stellar observations 



holds, one thing is clear: barring revolutionary advances 
in helio and asteroseismology, the information we have 
about solar and stellar magnetic activity is mainly lim- 
ited to each star's surface and atmosphere. Furthermore, 
it is the magnetic flux and energy that passes through 
the solar surface that shapes the structure and evolution 
of the solar corona and heliosphere and regulates space 
weather. Thus, understanding the link between dynamo 
action in the interior and flux emergence at the surface 
is a vital area of research for solar and stellar physics. 
While the solar dynamo operates on a wide variety of 
scales in both size and time, our simulations seek to make 
contact with elements of the large-scale magnetic behav- 
ior on time-sales of years to decades. 

1.2. Building Upon Theoretical Frameworks 

Cyclic dynamos are fundamentally three-dimensional, 
nonlinear, and chaotic. In spite of this difficulty, much of 
the groundwork for modern dynam o theory has been laid 
in analytic mean- field model s fe.g.. lParker|[l955t iMoffatl 
ll978tlKrause fe Radlerll 19801 ). The generation of toroidal 
field as differential rotation acts on a poloidal field, for 
example, can be well described using these models. The 
so-called i7-effect relies on shear from differential rota- 
tion in the convection zone or the tachocline at its base 
to stretch poloidal field into bands of toroidal field. The 
regeneration of poloidal field or the generation of oppo- 
site polarity poloidal field is parameterized in mean-field 
theory through the a-effect. A number of theoretical dy- 
namo models have been proposed, but as of yet no single 
numerical model has been able to capture all of the phys- 
ical m echanisms required (see review by iCharbonneaul 
[20H . 

To confront the complex nature of solar-like dynamo 
action, numerical models have been developed to explore 
aspects of vari ous d ynamo proc e sses. Pioneering work by 
I Oilman! (|1983f ) and lOlatzmaierl (| 19851) produced the first 
3D MHD simulation s of cyclic dynamo action in a rotat- 
ing spherical shell. iMiesch et all (120061 . I2008D have ex- 
plored the interplay of rotation, stratification, and mod- 
erately turbulent convection to produce strong differen- 
tial rotation in a hydrodynamic setting. When mag- 
netism is added in this context the resulting dynamo 
produces reversals of global polarity but is dominated 
by non-axisymme tric fields with little global organization 
(jBrun et al.ll2004 ). By adding an overshooting region of 
strong shear which mimics the solar tachocline, global- 
scale organization of the toroidal field was achieved by, 
but without reversals in global magnetic polarit y over 
about 30 si mulated years Browning et al.l (|2006[) . Re- 
cent work by iGhizaru et all (|2010D has shown large-scale 
organization of the toroidal field as well as magnetic ac- 
tivity cycles in a solar-like simulation; regular reversals 
of global magnetic polarity with a roughly 60 year pe- 
rio d for a comp l ete cy cle were achieved. Further work 
by iRacine et al.l (|2011l ) has interpreted these results in 
terms of mean-field dynamo theory. 

Additional insights into the solar dynamo have 
been realized through local simulations that do not 
model the full spherical geometry in order to achieve 
highe r resolution in a loca l domain (|Cline et al.l 
2003t iVasil fe Brummell I 2009D . The recent work of 



Guerrero fc Kapylal ( 2011 ) has shown that dynamo ac- 



tion in a domain with convection and a forced shear 
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Table 1 

Overview of Dynamo Cases 



Case 




Ra 


Ta 


Re 


Re' 


Rm 


Rm' 


Ho 


Roc 


V 


n 


Pm 


T E 


D3 


97 x 256 x 512 


3.28x10 s 


1.22X10 7 


173 


104 


86 


52 


0.374 


0.315 


13.2 


26.4 


0.5 


61.6 


D3a 


97 x 256 x 512 


5.84x10 s 


2.41 xlO 7 


244 


154 


122 


77 


0.447 


0.295 


9.40 


18.8 


0.5 


67.1 


D3b 


145 x 512 x 1024 


l.llxlO 6 


6.08X10 7 


343 


273 


171 


136 


0.566 


0.257 


5.92 


11.8 


0.5 


16.9 


D3-pml 


145 x 256 x 512 


2.98x10 s 


1.22X10 7 


149 


102 


149 


102 


0.372 


0.300 


13.2 


13.2 


1 


18.8 


D3-pm2 


145 x 512 x 1024 


3.08x10 s 


1.22X10 7 


145 


101 


291 


202 


0.370 


0.306 


13.2 


6.60 


2 


13.6 


S3 


145 x 512 x 1024 


7.68 x 10 s 


4.46 x 10 10 


8050 


5750 


4030 


2880 


0.581 


0.262 


0.218 


0.435 


0.5 


4.01 



Note. — Dynamo simulations at three times the solar rotation rate. All simulations have inner radius rb _ 
and outer radius of rt op = 6.72 X 10 10 cm, with L = (rt op — Hjot) = 1-72 x 10 10 cm the thickness of the spherical shell. Evaluated 
at mid-depth are the Rayleigh number Ra = (— 8p/8S)(dS /dr)gL 4 / puK, the Taylor number Ta = 4f2pL 4 /u 2 , the rms Reynolds 
number Re = v Yma L/v and fluctuating Reynolds number Re' = v' rlnB L/iJ, the magnetic Reynolds number Rm = v TulB L/r] and 
fluctuating magnetic Reynolds number Rm' = v' TIns L/rj, the Rossby number Ro = lj/2Qo, and the convective Rossby number 
Roc = (Ra/TaPr) 1 / 2 . Here the fluctuating velocity v' has the axisymmetric component removed: v' = V — (v), with angle 
brackets denoting an average in longitude. For all simulations, the Prandtl number Pr = u/k is 0.25 and the magnetic Prandtl 
number Pm = u/ij ranges between 0.5 and 4. The viscous and magnetic diffusivity, v and tj, are quoted at mid-depth (in units 
of 10 11 cm 2 s~ 1 ). The total evolution time Tg for each simulation is given in years. The values for case S3 with the dynamic 
Smagorinsky SGS model utilize the mean viscosity at mid-convection zone averaged on horizontal surfaces as well as in time. For 
case S3 using the dynamic Smagorinsky SGS model, the values quoted are based on the time-averaged rms viscosity, conductivity, 
and resistivity at mid-depth, noting that these diffusion coefficients have near hundred-fold spatial variations. 



layer can produce strong dynamo action and yield buoy- 
ant magnetic structures. Another approach using helical 
forcing in a Cartesian domain has shown that large-scale 
magnetic structures which undergo regular reversals in 
polarity can be achieved even without conv ection or ro- 
tatio n, hinting at the key role of helicity (jMitra et al.1 
f2010h . 

He re we expand upon the work of iBrown et al.l (|2010l 
120 111 ) in exploring global 3D simulations of dynamo ac- 
tion in sun-like stars rotating faster than the current solar 
rotation rate. We report on a suite of simulations at 3J1q 
which explore the transition from dynamos with persis- 
tent toroidal fields to cyclic dynamos by decreasing the 
level of explicit diffusion in the simulations. Although 
the simulations reported here are ostensibly at a rotation 
rate of 3^0, it is important to realize that the dynamics 
we describe may not be restricted to young stars. With 
regard to the generation of magnetic wreaths, the most 
salient nondimensional parameter of the physical system 
is the Rossby number, Ro= cj/(2f2o)> which is small in 
both the Sun and these models. In this series of models, 
we begin to address the feasibility of the wreath-building 
dynamo mechanisms at higher levels of turbulence and 
study the global-scale reversals and cycles of magnetic 
activity achieved, which will be the focus of §3-6. We 
also begin to explore the subtle but crucial link between 
wreath-building convective dynamos and flux emergence, 
which will be the focus of §7. 

2. DYNAMOS AT 3Q Q 

We study convection and dynamo action in the deep in- 
terior of solar-like s tars using the anelastic spherical har - 
monic (ASH) code (IClune et al.lll999t IBrun e"t~aT1[200l . 
Our simulation approa ch is briefly describe d here, but is 
more fully explained in IBrown et al.l (|2010l ). ASH solves 
the anelastic MHD equations in a rotating spherical shell 
with a background stratification taken from a ID model 
of solar structure. We focus on simulating the bulk of 
the solar convection zone from 0.72i? Q to O.97i?0 (i?© 
is solar radius) with a density contrast of about 25. We 



do not model the near-surface layers of the sun, for we 
are limited by the anelastic approximation to subsonic 
flows. Additionally we cannot resolve the small-scales 
of motion needed to simulate granular and supergranu- 
lar scales. We also do not include the stably-stratified 
radiative zone or the tachocline in these simulations, al- 
though simulations includi ng those compone nts are an 
active area of research (see IBrun enu1l2CLTTI) . We have 
done some preliminary work in adding a tachocline to 
these simulations and found that it does not drastically 
change the dynamo action in the bulk of the convective 
layer. The effects of a tachocline will will explored further 
in a future paper. Our results tend to support the recent 
studies with mean-held dynamo models, which suggests 
that the differential rotation of the convection zone may 
play a greater role in the gener ation of toroidal magnetic 
field than the tachocline (e. g., iDikpati fc Gilmanl 120061 : 
iMuhoz Jaramillo et al.ll2009f ). We use impenetrable and 
stress-free boundary conditions on both the top and bot- 
tom of the domain. We impose the entropy gradient at 
the top and bottom of the domain for the thermal bound- 
ary conditions. For the magnetic fields we use a perfect 
conductor condition on the bottom boundary and match 
to an external potential field on the top boundary. These 
conditio ns and our evolution equations are described in 
detail in IBrown et al.l (|2010D . 

ASH is a large-eddy simulation which employs a 
subgrid-scale model to account for the effects of un- 
resolved scales of motion. The standard subgrid-scale 
(SGS) model in ASH simulations uses enhanced values 
of viscosity, thermal diffusivity, and magnetic resistivity 
relative to those expected from atomic values in order 
to represent additional mixing due to unresolved turbu- 
lent motions. In this enhanced eddy SGS model, viscos- 
ity v, thermal diffusivity n, and magnetic resistivity 77 
all scale as p -1 ^ 2 , where p is the spherically-symmetric 
background density of the simulation. This prescription, 
along with constant Prandtl and magnetic Prandtl num- 
bers throughout the domain, follows that of lBrown et al.l 




Figure 1. Magnetic wreaths, (a)-(c) Shown in global Mollweide view (equator at middle, poles at top and bottom) is the radial velocity 
of the convection at 0.95Rq in cases D3, D3a, and D3b, respectively, (d)-(f) Also in Mollweide few, longitudinal magnetic field at 
mid-convection zone at times ti indicated in Figure [3] (g)-(i) Shown at the same times for each case is a 3D field line rendering of the 
magnetic wreaths near the equator. In both types of display for the magnetic field, color gives the polarity and amplitude of the longitudinal 
field (red positive, blue negative). Times shown correspond to t\ for each case in Figure [3] 



(l2010l[20Tll) . All cases presented in this paper use Pr = 
v/k = 0.25, but variable Pm (see Table [TJ. 

In addition, we have also implemented a more complex 
SGS tre atment, the dyn a mic S magorinsky model devel- 
oped by IGermano et all (|1991l ). By using the dynamic 
Smagorinsky model in ASH simulations we are able to 
reduce the mean diffusion at mid-convection zone by a 
factor of 50 without an increase in resolution. Our imple- 
mentation of the dynamic Smagorinsky model is summa- 
rized in Appendix A. This SGS trea tment is only used i n 
case S3, which was first presented in lNelson et al.l (| 2 1 If ) . 

Table Q] presents the computational resolution, rele- 
vant non-dimensional parameters, diffusion coefficients, 
and total evolution time for each of the six cases we will 
discuss here. We have explored two main branches in pa- 
rameter space. The first branch includes cases D3, D3a, 
and D3b, where viscosity v, thermal diffusivity n, and 
magnetic resistivity 77 have all been dropped together, 
thus keeping a constant magnetic Prandtl number. The 
second branch includes cases D3, D3-pml, and D3-pm2, 
where v and n are held constant and only 77 is decreased, 
resulting in increasing magnetic Prandtl numbers. We 
will refer to the two branches as the constant Pm and 
increasing Pm branches, respectively. The constant Pm 
branch was found to be more compelling, as cases D3a 
and D3b generally produced strong magnetic wreaths 
that were anti-symmetric about the equator, whereas 
the high Pm branch produced a wider variety of sym- 
metric and anti-symmetric toroidal field states and was 
therefore less amenable to study. Such behavior is not 
unexpected as dynamos with higher magnetic Prandtl 
number tend to promote small-scale dynamo action. Wc 



will generally focus on the constant Pm branch of sim- 
ulations while referencing the increasing Pm branch to 
provide additional insight. 

Case D3 was initiated from a well developed hydro- 
dynamical simulation that was seeded with a small ran- 
dom magnetic field. Each subsequent case along both 
branches was started from the preceding case. Thus both 
cases D3a and D3-pml were started using case D3 as ini- 
tial conditions, case D3b was started using case D3a, and 
so on. Wc have re-started case D3a from a random seed 
field to verify that it settles into a similar region of so- 
lution space as the version continued from case D3 and 
found no strong differences in the time-averaged behavior 
over several thousand days. 

3. MAGNETIC WREATHS 

The dominant magnetic structures built by each of 
these simulations are the low latitude bands of predom- 
inately toroidal field, which we term wreaths. These 
wreaths are generally anti-symmetric about the equa- 
tor, though symmetric states are observed along with 
states where one hemisphere displays a wreath while the 
other does not. These irregular states are most com- 
mon along the increasing Pm branch of our simulations. 
The wreaths in ca se D3 are discussed extensively by 
I Brown et~a l. (2010) and additional wreaths are analyzed 
at som ewhat faster rotation rate (5f2©) by IBrown et al.l 

pal . 

3.1. Magnetic Topology 

Figure Q] shows snapshots of the turbulent convection 
and the wreaths for cases D3, D3a, and D3b at mid- 



Magnetic Wreaths and Cycles 



■5 



Table 2 

Volume- Averaged Energy Densities and Differential Rotation Rates 



Case 


Total ME 


TME 


PME 


FME 


Total KE 


DRKE 


MCKE 


FKE 






D3 


0.68 (9%) 


0.29 (43%) 


0.029 (4%) 


0.36 (53%) 


6.67 (91%) 


4.35 (65%) 


0.010 (0.1%) 


2.31(35%) 


112 


192 


D3a 
D3b 


0.88 (12%) 
0.82 (13%) 


0.32 (36%) 
0.10 (12%) 


0.030 (3%) 
0.011 (1%) 


0.52 (59%) 
0.70 (85%) 


6.41 (88%) 

5.42 (87%) 


3.71 (58%) 
2.45 (45%) 


0.011 (0.2%) 
0.012 (0.2%) 


2.68 (42%) 
2.96 (55%) 


101 

95 


163 
131 


D3-pml 
D3-pm2 


1.04 (18%) 
1.17 (21%) 


0.26 (25%) 
0.15 (13%) 


0.033 (3%) 
0.028 (2%) 


0.75 (72%) 
0.99 (85%) 


4.87 (82%) 
4.34 (75%) 


2.63 (54%) 
2.29 (53%) 


0.010 (0.2%) 
0.009 (0.2%) 


2.23 (46%) 
2.04 (47%) 


87 
74 


139 
121 


S3 


0.83 (13%) 


0.072 (9%) 


0.0060 (0.7%) 


0.75 (91%) 


5.50 (87%) 


2.32 (42%) 


0.013 (0.2%) 


3.17 (58%) 


95 


133 



Note. — Volume-averaged magnetic and kinetic energies for dynamo simulations at three times the solar rotation rate, as well as the 
magnitude of the differential rotation contrast in radius at the equator Af2 r and the average contrast at the top of the simulated domain 
between the equator and ±60° latitude AQg. Shown in units of 10 6 erg cm -3 are the total magnetic energy (Total ME), axisymmetric 
toroidal magnetic energy (TME), axisymmetric poloidal magnetic energy (PME), fluctuating magnetic energy (EME), total kinetic energy 
(Total KE), differential rotation kinetic energy (DRKE), meridional circulation kinetic energy (MCKE), and fluctuating kinetic energy 
(FKE). The percentage of the total energy is shown for total magnetic energy (Total ME) and total kinetic energy (Total KE). The 
percentage of the total magnetic or kinetic energy for each component is shown in parentheses. Values for differential rotation rates are 
in units of nHz (3^0 = 1240 nHz). Values are averaged in time over long intervals. 



convection zone in global Mollwcidc view as well as at 
low latitudes in a 3D volume rendering of magnetic field 
lines colored by B^. In all three cases strong longitude- 
averaged fields arc present at low latitudes, however the 
nature of the wreaths change from case D3 where axisym- 
metric fields dominate to case D3b where a significant 
axisymmetric field component is present but not domi- 
nant. In case D3b the morphology has changed such that 
the wreaths are confined in longitudinal extent. Figure 
[1] shows a typical field configuration, but the wreaths are 
observed at various times to extend over as little at 45 de- 
grees and as much as 270 degrees in longitude. All three 
cases show extensive connectivity between the wreaths 
and the surrounding domain where magnetic fields are 
moderate in strength but far less coherent. The wreaths 
are strongly modulated by the convective flows, produc- 
ing a ragged appearance that is particularly noticeable 
in case D3b but present in all three cases. In the more 
turbulent cases there are also significant small-scale mag- 
netic fields at moderate to high latitudes, and occasional 
locally-generated wreath-like structures near the poles 
which persist for less than about 100 days at a time. 

The shift from structures dominated by axisymmet- 
ric fields in case D3 to the patchy wreaths in case D3b 
is illustrated by the changes in the time- and volume- 
averaged energy densities shown in Table [5J Between 
cases D3 and D3b there is a roughly 30% increase in 
the total magnetic energy of the simulation, though 
the energy in both the axisymmetric toroidal (TME) 
and poloidal (PME) fields decreases by roughly a fac- 
tor of three. The doubling of the energy in the non- 
axisymmetric magnetic fields more than compensates for 
the decrease in mean fields. When compared with the 
kinetic energy densities, the changes in the magnetic en- 
ergies becomes even more striking. Viscous, thermal, 
and magnetic diffusion in case D3b are all reduced by 
the same factor relative to case D3. However the to- 
tal kinetic energy in case D3b dropped by 19%. The 
non-axisymmctric kinetic energy (FKE) rose only mod- 
erately compared to the decrease in differential rotation 
kinetic energy (DRKE). The high magnetic Prandtl num- 
ber cases also show a tendency towards larger total and 
fluctuating magnetic energies, as well as reduced axisym- 
metric toroidal magnetic energy as the magnetic diffusion 




Figure 2. Probability distribution functions for unsigned B^, at 
mid-convection zone for cases D3 (purple), D3a (green), D3b (red), 
and S3 (blue) showing the surface area covered by fields of a given 
magnitude. Distributions are averaged over about 300 days when 
fields are strong and as steady as possible. Dashed vertical lines 
show the field-strength at which cquipartition is achieved with the 
maximum fluctuating kinetic energy (FKE) at mid-convection zone 
for each case. Case D3b shows a deficit of field in the 10 kG range, 
but an excess of surface area covered by extremely strong fields 
above 25 kG range, as well as higher peak field strengths. Case S3 
shows significantly greater regions of fields in excess of 20 kG than 
all other cases. 



is reduced. 

It is illustrative to compare cases D3b and D3-pmf, 
as they have roughly equal levels of magnetic diffusion, 
with case D3b having comparatively lower levels of vis- 
cosity and thermal diffusion. The largest differences are 
in the axisymmetric magnetic energies which are both 
about three times greater in case D3-pml than in case 
D3b. This may be due to the more laminar flow in case 
D3-pmI, which would tend to create fewer sharp gradi- 
ents in the large-scale magnetic structures and thus lower 
the effective dissipation in case D3-pmI compared to case 
D3b, even though the diffusion coefficients in the induc- 
tion equation are nearly the same. Case D3-pml also 
show significantly less differential rotation contrast both 
in radius and latitude compared to case D3b, pointing to 
the key role of magnetic torques in decreasing differential 
rotation, which will be discussed further in 55. 
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3.2. Non-axisymmetric Fields 

Our discussion of the magnetic wreaths to this point 
has focused on the axisymmetric fields, which are pro- 
gressively weaker in moving from case D3 to case D3b. 
While the axisymmetric fields weaken with increased tur- 
bulence, very strong fields become more common when 
measured by the fraction of the domain they occupy. Fig- 
ure [2] shows the probability distribution function for 
at mid-convection zone in cases D3, D3a, D3b, and S3. 
While case D3b has a deficit of fields around 10 kG com- 
pared to case D3a, there is a clear excess of fields above 20 
kG. Interestingly the distribution for case D3b is greater 
than that for case D3 for all but the smallest bin, indicat- 
ing that while case D3 may have stronger axisymmetric 
fields in the low latitude wreaths, case D3b compensates 
by having higher amplitude fluctuating fields throughout 
the domain. The peak field strength at mid-convection 
zone is 32 kG in case D3, 36 kG in case D3a, and 38 kG 
in case D3b. Near the base of the convection zone case 
D3b exhibits even stronger fields of up to 44 kG. Case S3 
posses magnetic fields of up to 45 kG at mid-convection 
zone and 52 kG near the base of the convective layer. For 
all four cases fields are seen well in excess of equipartition 
energies with the maximum fluctuating kinetic energy of 
the flows. This is a clear indication of turbulent inter- 
mittency in the magnetic fields. 

A statistical measure of turbulent intermittency is the 
time-averaged excess kurtosis given by 



(a; 



kurt{B } 



I-M-BA f [B'A dB' 4 



3, (1) 



where f(B'i) is the probability distribution function (see 
lPopd[2000l ). For reference a Gaussian distribution would 
have an excess kurtosis of 0. The level of turbulent in- 
termittency is measured by how leptokurtic the distribu- 
tion is found to be, with large values corresponding to 
increased intermittency. For case D3 kurtjB^} = 9.6, 
while for case D3a kurtjiJ^} = 10.5, and for case D3b 
kurtj-B^} = 12.1. Leptokurtic distributions are likely 
to experience strong coherent structures, such as the 
strong regions of coherent toroidal field in these simu- 
lations. At even lower levels of diffusion than can be 
realized with the enhanced eddy SGS model, the strong- 
field regions become sufficiently buoyant and coherent so 
as to form buoyant magnetic loops as realized in case S3 
(jNelson et al.ll2MTl ). for which kurt{B } = 15.6. Highly 
leptokurtic distributions like these indicate that extreme 
events are enhanced relative to a Gaussian distribution, 
and the trend towards increasing kurtosis as simulations 
become more turbulent points to turbulent amplifica- 
tion of magnetic fields. As we will discuss further in 
§7, this provides an alternate pathway to produce re- 
gions within the larger wreaths which can be amplified 
through turbulent intermittency to produce coherent re- 
gions of strong magnetic field, which can then become 
buoyant. We term this the turbulence-enabled magnetic 
buoyancy paradigm. 

4. CYCLIC REVERSALS ACHIEVED BY REDUCING 
DIFFUSION 
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Figure 3. Time-latitude plots of longitude averaged toroidal mag- 
netic field (B^) at O.79i?0 for (a) case D3 over about 56 years, (b) 
case D3a over the same amount of time, and (c) case D3b over 
ab out 13 years. Dotted lines show times referenced in Figures [T] 
[TJ 1111 1121 and 1151 Dashed lines on (b) indicate the time period 
covered by (c). Case D3b was started from case D3a at <2 (dotted 
line). The evolution of case D3b is limited by the increased com- 
putational cost of the higher resolution required for computational 
stability. 



In addition to building strong magnetic wreaths, cases 
D3a and D3b exhibit cyclic reversals of global magnetic 
polarity. As is believed to occur in the sun, the gen- 
eral pattern of the cycles is that the toroidal fields peak 
at roughly the time when the poloidal field is reversing 
sign, and the poloidal fields peak in amplitude when the 
toroidal fields are reversing sign. There are also a number 
of variations on this pattern, where one hemisphere may 
develop considerably stronger fields than the other or 
where both hemispheres have the same sense of toroidal 
field, pointing to large contributions at these times from 
quadripolar poloidal fields. Cases D3-pml and D3-pm2 
also display strong variations in the strength and topol- 
ogy of their axisymmetric fields. However the irregular- 
ities are more pronounced for these cases over the time 
simulated. 

4.1. Reversals in Global Magnetic Polarity 

Figurc[3]shows the temporal evolution of the longitude- 
averaged toroidal field {B^} at mid-convection zone over 
the history of cases D3 (Figure Ha)), D3a (Figure fflp)), 
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Figure 4. Reversal in magnetic polarity of the toroidal wreaths in 
case D3b shown in in Mollweide projection at mid-convection 
zone on left, and in (B^) in longitudinal average over latitude and 
radius on right. Color indicates strength of toroidal magnetic field 
with the color table saturating at ±7 kG for the Mollweide im- 
ages and ±3 kG for the longitudinal averages. Five snapshots 
corresponding to t2 through t& from Figure (3{c) arc shown each 
separated by roughly 120 days. 



and D3b (Figure EJc)). In case D3 we see persistent 
wreaths centered at about 20° above and below the equa- 
tor. These wreaths persist for about 68 years or as long as 
we have run the simulation. The polarity of the wreaths 
is constant in time, though variations on roughly 6 year 
time scales can be seen in both the amplitude of the low 
latitude wreaths as well as the propagation of field to 
higher latitudes. The behavior of this case is discussed 




Time (years) 

■10 KG^BZZ^MIO kG 

Figure 5. (a) Time evolution in case D3-pml of the axisymmet- 
ric toroidal magnetic field at 0.79 Rq over roughly 15 years of 
simulated time. Strong variability of the mean fields is seen in 
both hemispheres, (b, c) Companion snapshots of B^ at 0.84 Rq 
showing the spatial variability and non-axisymmctric nature of the 
wreaths. Successive snapshot times are indicated by dashed lines 
in (a). 



in detail in iBrown et al.1 (|2010f ). Figure El[b) shows case 
D3a over a comparable length of time as in the first panel. 
Case D3a undergoes reversals in global magnetic polarity 
as well as three significant irregular states. Additionally 
there arc modulations in the amplitude of the wreaths 
and poleward movements of field on roughly 3 year time 
scales. These variations are not always synchronized be- 
tween the two hemispheres, and neither are the reversals, 
indicating that the poloidal field can have a complicated 
structure. 

Figure [3{c) shows the temporal evolution of {B$) at 
mid-convection zone for case D3b over about 13 years, 
with indications of cycles of magnetic activity and rever- 
sals of global polarity. We have simulated 10 reversals as 
measured by the time-smoothed antisymmetric compo- 
nent of (B$) changing sign. The time between reversals 
ranges from 0.6 to 1.9 years and, as in case D3a, the 
two hemispheres are not always synchronized. There are 
several times when one hemisphere shows significantly 
stronger fields than the other or when both hemispheres 
have the same sense of fields. This is partly due to the 
averaging procedure used to create these figures and the 
fact that we are only looking at a single depth. Analysis 
of the full 3D data shows that there is almost always a 
wreath-like structure in each hemisphere. 

Figure |4] shows a sequence of snapshots of B§ at mid- 
convection zone in case D3b over a full spherical shell 
and of (Be/,) over the domain before, during, and after a 
reversal in the polarity of the wreaths. Each snapshot 
is roughly 120 days after the previous snapshot. The 
wreaths start appearing as strong mean-field structures 
in the longitudinal average, but the non-averaged cut 
at mid-convection zone shows that there is significant 
longitudinal variation in the wreaths, with the northern 
wreath covering roughly 120° and the southern wreath 
covering 180° in longitude. There is also substantial 
evidence for interactions between the wreaths near the 
center of the image in Figure SJa). As time progresses, 
the axisymmctric fields weaken as the non-axisymmctric 
components begin to dominate. Small patches of strong 
field persist, but they are largely washed out in the lon- 
gitudinal averages. After about 480 days (Figure H]i) 
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Figure 6. (a) Time evolution in case D3-pm2 of the axisymmctric 
toroidal magnetic field at 0.79 Rq over roughly f3 years of sim- 
ulated time. Strong variability of the mean fields is seen in both 
hemispheres, along with irregular reversals in polarity, at times 
in only one hemisphere and at other times globally, (b, c) Com- 
panion snapshots of at 0.84 Rq showing the spatial variability 
and non-axisymmetric nature of the wreaths. Successive snapshot 
times are indicated by dashed lines in (a). 



strong patches of opposite polarity field begin to appear 
and by the final frame (Figure 2^) the strong mean fields 
have been reestablished in the opposite hemispheres from 
the initial configuration. 

4.2. Variability at Higher Magnetic Prandtl Number 

The simulations on the increasing Pm branch also show 
increased temporal variability relative to case D3. There 
is also evidence for a change in the nature of the dynamo 
action in these simulations. Figure [5] shows the evolu- 
tion of (B^) at 0.79 Rq over the history of case D3-pml, 
along with snapshots of the toroidal magnetic field at 
mid-convection zone at three representative times. This 
case has selected a configuration of toroidal field that is 
largely symmetric about the equator and of essentially 
the same polarity at most times. Some periods of pos- 
itive polarity field are seen, though the dominant field 
in both hemispheres is clearly of negative polarity. Un- 
like cases D3a and D3b, case D3-pml does not undergo 
a true global reversal of magnetic polarity. It does, how- 
ever, exhibit strong temporal variability in the wreaths 
seen in both hemispheres to an extent not seen in cases 
D3 or D3a. 

Figure |6] shows a similar view of case D3-pm2 over its 
temporal evolution. Again it tends to avoid the anti- 
symmetric states characteristic of cases D3, D3a, and, to 
a lesser extent, D3b. This case, however, docs exhibit 
clear reversals of global magnetic polarity. Interestingly, 
these reversals do not appear to occur at regular inter- 
vals, and often one hemisphere can reverse without a no- 
ticeable change in the other hemisphere. As an example, 
the southern hemisphere maintains a positive polarity 
wreath between about t = 5.5 years and t = 10.5 years 
while the northern hemisphere exhibits four reversals in 
that same time interval. 

The preference for irregular polarity states in B^ along 
the increasing Pm branch is clearly related to the de- 
creased level of magnetic diffusion, though it may also 
be indicative of a shift in behavior due to the transition 
from small to large magnetic Prandtl number. In cases 
D3, D3a, and D3b magnetic diffusion occurs on scales 
larger than those related to the diffusion of momentum. 



This tends to promote the concentration of magnetic en- 
ergy at large scales. For high Pm dynamos, the resistive 
scale is smaller than the viscous scale, which tends to 
prom ote the growth of magne tic energy at small scales 
(e.g.. ISchekochihin et al.ll2004T ). There is still consider- 
able large-scale organization of magnetic field by the dif- 
ferential rotation, but the increasing Pm branch exhibits 
less ordered behavior than the constant Pm branch of 
simulations. 

When examining the relative importance of decreased 
magnetic diffusion and increased levels of turbulence, it 
is perhaps most instructive to compare cases D3b and 
D3-pm2. Table [5] shows that the division of magnetic 
energies between the axisymmetric toroidal, axisymmct- 
ric poloidal, and fluctuating magnetic fields are roughly 
equivalent in the two cases, although case D3-pm2 has 
more magnetic energy overall. The kinetic energies in 
case D3-pm2 are, however, more similar to case D3 than 
case D3b with the exception of decreased differential ro- 
tation kinetic energy due to enhanced Lorentz force feed- 
backs. This suggests that the onset of reversals is driven 
primarily by decreasing magnetic diffusion rather than 
by some subtle shift in the velocity fields or correlations 
between magnetic fields and velocities on small scales. 

5. MAINTAINING ROTATIONAL SHEAR 

A crucial component in the construction of magnetic 
wreaths is the strong latitudinal and radial shear from 
the differential rotation. The ^-effect has previously 
been shown to be the primary production mechanism for 
the magnet ic wreaths in cases D3 and D5 (jBrown et al.1 
I2010L 120 111 ), and it plays a key role in these simula- 
tions as well. Thus the angular momentum transport 
required to maintain the differential rotation is an im- 
portant physical process in these dy namo models . In th e 
hydrodynamic models explored by iBrown et al.l (|2008| ). 
angular momentum transport in simulations at 3f2© was 
shown to be a balance between Reynolds stress support- 
ing solar-like differential rotation with the meridional cir- 
culation and viscous diffusion tending to dissipate gradi- 
ents in the rotation profile. With the addition of mag- 
netic fields, Maxwell stress and mean magnetic torques 
can also transport angular momentum, changing the bal- 
ance supporting the strong differential rotation achieved 
in the hydrodynamic cases. Even in cases witho ut mag- 
netic cycles such as case D3. IBrown et~aT1 (|2010[ ) showed 
that there are significant feedbacks on the differential 
rotation profile due to variations in the strength of the 
magnetic fields over time. It is thus useful to examine 
not only the steady state balance of angular momentum 
transport over long time averages covering many mag- 
netic cycles, but also to look at the temporal variability 
of those balances. 

In order to explore the transport of angular momen- 
tum, let us examine the physical mechanisms which come 
into play. The balance of specific angular momentum 
along the rotation axis is determined by taking the prod- 
uct of the cylindrical radius A = r sin 9 and the longitu- 
dinal component of the longitude-averaged momentum 
equation, which can be expressed as 



at 



V ■ T. 



(2) 



We decompose the flux vector of mean angular momen- 
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Figure 7. Differential rotation and the terms contributing to the accompanying redistribution of angular momentum in case D3b averaged 
over many magnetic cycles, (a) Angular velocity Q profile with radius and latitude, accompanied in turn by profiles of specific angular 
momentum flux in cylindrical radius (A) by Reynolds stress (RS), meridional circulation (M C), vis cous diffusion (VD), Maxwell stress 
(MS), and mean magnetic torques (MT), respectively. Terms are defined in detail in Equations IB II and lB2l They are here averaged in time 
and longitude, and given in units of 10 15 g s — 2 . (b) Scalar plot of 2-integrated fluxes of angular momentum with cylindrical radius A in 
units of 10 38 g cm 2 s — 2 . Reynolds stress (RS, red) balance Maxwell stress (MS, blue), with viscous diffusion (VD, green) and magnetic 
torques (MT, purple) playing less of a role. Contribution from the meridional circulation (MC, brown) are small. The sum of all five terms 
are also plotted (black dashed line). 



Table 3 

Production and Dissipation of Differential 
Rotation Kinetic Energy 



Case 


Lrs 


Lmc 




Lms 


Lmt 


D3 


4.26 


-0.020 


-2.45 


-1.36 


-0.68 


D3a 
D3b 


3.18 
2.59 


-0.032 
-0.003 


-1.11 
-0.64 


-1.36 
-1.94 


-0.70 
-0.19 


D3-pml 
D3-pm2 


3.64 
3.32 


-0.014 
-0.024 


-1.68 
-1.61 


-1.87 
-1.95 


-0.35 
-0.31 



Note. — Total production and dissipation 
of kinetic energy in the axisymmctric differen- 
tial rotation profile over the entire simulated 
volume and averaged in time. Values for en- 
ergy production rates are in units of 10 32 erg 
s — 1 . Production and dissipation terms are split 
following Equation ((3j into contributions from 
viscous dissipation, Reynolds stress, meridional 
circulations, Maxwell stress, and mean magnetic 
torques, respectively. Production terms are de- 
fined in Appendix B. 




D3-Pm2 D3-Pml D3 D3a D3b 



Figure 8. Companion to Table \3\ showing the balance of time- 
averaged generation terms for the kinetic energy in the differential 
rotation profiles for each case indicated. In all cases the differential 
rotation is maintained by a balance between the Reynolds stress 
and a combination of viscous diffusion and fluctuating and mean 
magnetic torques. The contribution from meridional circulations 
are not shown due to their small magnitude. 



turn T into radial and latitudinal components followin 
prior convention (lElliott et all 120001: IBrun et all 1200 
iBrown et al.l 120 111 ) . We also decompose the flux vec 
tor into cylindrical coordinates along cylindrical radius 
(A) and along the rotation axis (z), which in many ways 
is advantageous for displaying these quantities. A de- 
tailed description of this decomposition is given in Ap- 
pendix B. The cylindrical flux of angular momentum is 
shown for case D3b over a long time average in Figure [3 
The differential rotation is again clearly maintained by 
the Reynolds stress (RS), however here the terms oppos- 
ing the differential rotation have changed compared to 
similar hydrodynamic cases. In case D3b the Maxwell 
stress (MS) is the largest term opposing the Reynolds 
stress with viscous diffusion (VD) and the mean mag- 
netic torques (MT) each playing a small role, while con- 
tribution of the meridional circulation (MC) is almost 
insignificant. 

We can write the evolution of the total energy of the 



differential rotation £dr as 
<9£dr 



dt 



= L 



VD 



L 



RS 



L 



MC 



L 



MT, 



(3) 



where the terms on the right-hand side represent the 
sources and sinks of kinetic energy in the differential ro- 
tation due to, respectively, viscous diffusion, Reynolds 
stress, meridional circulations, Maxwell stress, and mean 
magnetic torques. Appendix B provides a derivation of 
Equation ([3]) and an expansion of the sources and sinks. 

Using this decomposition, we can examine the bal- 
ance of production and dissipation of £dr averaged over 
long time intervals in each simulation. The balances are 
represented in Table [3] and Figure [8] For the increas- 
ing Pm branch the Reynolds stress change only slightly 
while the mean magnetic torques and viscous diffusion 
are systematically replaced by the Maxwell stress. Sim- 
ilar trends are observed in the constant Pm branch of 
cases, though here the magnitude of the Reynolds stress 
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Figure 9. Temporal variability o f diff erential rotation in case D3b 
over the same interval as in Figure ITB1 (a) Longitudinally-averaged 
rotation rate at mid-convection zone as a function of time and 
latitude, (b) Temporal variations are accentuated by subtracting 
the time-averaged Q at each latitude. Bands of faster rotating fluid 
move poleward on about the cycle period, (c) Rotation contrasts 
in radius at the equator A r (red, solid) and in latitude between the 
equator and ±60° in the upper convection zone Ag (red, dashed). 
The volume-averaged toroidal field strength is also shown (blue, 
solid), with a phase-lag between peaks in magnetic field strength 
and decreases in differential rotation. Dotted lines indicate times 
t2 through tg from Figure [4] 

and viscous diffusion terms decrease more dramatically. 
This shift from unresolved dissipation in the form of 
SGS viscosity to resolved, small-scale torques from the 
Maxwell stress indicates that the balances which main- 
tain strong differential rotation can persist in less dif- 
fusive regimes, assuming that magnetic energies remain 
significantly smaller than kinetic energies. 

Turning to the temporal variability in these balances, 
we find that for case D3b the departures from the values 
presented in Tableland Figurc[8]are about 10% for Lms 
and Lmt when averaged over about 10 days, whereas 
those in Lrs, Lmc and Lvd arc about 1%. This leads 
to decreases in the differential rotation when magnetic 
fields are strong, such as near the peak of the magnetic 
activity cycles. Conversely, we observe modest increases 
in the differential rotation when magnetic fields are weak, 
such as during reversals of magnetic polarity. 

Figure [5] shows the differential rotation in case D3b 
over several magnetic reversals. The differential rota- 
tion profile at mid-convection zone in Figure EJa) is per- 
sistent, though there are small systematic variations in 
f2 revealed in Figure EJb) during each magnetic cycle. 
Figure [9jc) shows the differential rotation contrast in 
both radius and latitude over time as well as the volume- 
averaged toroidal magnetic field strength, indicating that 
the modest variations in the differential rotation are re- 
lated to those in the magnetic field. 

6. PRODUCTION OF MAGNETIC FIELD 

The transition from persistent wreaths in case D3 to 
cyclic wreaths and global polarity reversals in case D3b 
indicates that by reducing the levels of diffusion in these 
simulations we have fundamentally altered the balance 



of terms in the magnetic induction equation. The details 
of the reversal mechanism are likely to be very subtle in 
these highly nonlinear systems. In order to better under- 
stand the reversal mechanism, we explore the nature of 
the balances in the production and dissipation of toroidal 
and poloidal magnetic fields and provide indications of 
where and why changes in those balances are occurring, 
as well as some hints as to the nature of the reversal 
mechanism. 

6.1. Generation of Toroidal Magnetic Energy 

In iBrown et al.1 (|2010T l a detailed analysis of the bal- 
ance of toroidal component of the axisymmctric induc- 
tion equation was presented. We write the toroidal com- 
ponent of the induction equation as 

f V x (v x b)} -Vx(,,Vx (B^j) • (4) 



dt 

Using vector identities, the first term on the right-hand 
side can be written as the sum of shearing terms, ad- 
vection terms, and a compression term; additionally all 
of these terms can be decomposed into mean and fluc- 
tuati ng components (for a full derivation, see Appendix 
A in IBrown et al.l I2010D . That work also showed that 
the wreaths in case D3 are primarily generated by the 
J7-effect and dissipated by a combination of small-scale 
advection, shear, and diffusion. 

Here we perform a similar analysis, but instead of ex- 
amining the generation of {B^), we choose to examine 
the generation of the volume-integrated energy of the ax- 
isymmctric toroidal fields over the entire computational 
domain V, given by 



£tme 



dV. 



(5) 



We can construct an evolution equation for £tme by mul- 
tiplying Equation ((4]) by (B^). The result can be written 



d& 



TME 



dt 



Cms + Gfs + Gma + Gfa + Gac + Grd, (6) 



where the six terms on the right-hand side represent, 
from left to right, the shearing of axisymmctric mag- 
netic fields by mean flows associated with the ^-effect 
(Cms), the average of fluctuating flows shearing fluctua- 
tion fields (Gfs)j the advection of mean fields by mean 
flows (Gma), the average of fluctuating flows advert- 
ing fluctuating fields (Gfa), the anclastic compression 
of fields (Gac), and the resistive diffusion of mean fields 
(Gr,d). Unlike in previous analyses which looked at the 
generation of magnetic field vectors, here we are con- 
cerned with scalar quantities. The terms on the right- 
hand side of Equation © are computed as 



Gms = / {B4) [((B) -v) {v)] dV 



Gfs : 
Gma 
Gfa = 



(B4,) 



B+) {(v)-V)(B) 



(B^) ((v'-V)B' 



dV, 
dV, 
dV, 



(7) 
(8) 
(9) 
(10) 
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Figure 10. Volume-integrated production terms of magnetic en- 
ergy in the mean toroidal fields from Equation (JBJl for (a) case 
D3 and (b) case D3b. We have combined the mean shear and 
mean advection terms (blue line) and the fluctuating shear and 
fluctuating advection (purple line). In both simulations, energy is 
produced primarily by the shearing of mean fields by mean flows. 
Also in both cases compression of fields (green line) plays a very 
small role. In case D3 diffusion (red line) and the advection and 
shear of fluctuating flows on fluctuating fields destroy energy, with 
diffusion generally a factor of 2.5 larger. In case D3b, however, 
the dissipation of energy by fluctuating advection and shear is 2.2 
times greater on average than diffusion. Thus in case D3b resolved 
turbulence is the primary mechanism for dissipating the magnetic 
wreaths. 




Grd = - / (B )V x (77V x dV. (12) 

Jv 

For consistency, angle brackets denote longitude aver- 
ages. 

Figure [TU] shows the temporal evolution of the inte- 
grated energy of the axisymmetric toroidal field for cases 
D3 and D3b and the behavior of the production terms 
governing the variation of £tme- We have chosen to com- 
bine the contributions of the mean shear and advection 
terms and the fluctuating shear and advection terms for 
ease of viewing. The mean advection term Gma is gen- 
erally positive and always much smaller than the mean 
shear term Gms- The fluctuating shear and advection 
terms are both generally negative, of approximately the 
same magnitude, and tend to vary in phase with each 
other. 

Let us first look at the average levels of each term plot- 
ted in Figure [10] to get a sense for the basic balance of 
terms. The production of £tme is dominated by the 



mean shear term which is large and always positive in 
both case D3 and D3b. The compression term in both 
cases is roughly an order of magnitude smaller but is 
again always positive due to the asymmetry in upflows 
and downflows in compressible convection, which gives 
preference to downward pumping of magnetic field caus- 
ing an increase in magnetic energy due to compression. 
The production of magnetic energy is opposed by the 
resistive diffusion, fluctuating advection, and fluctuating 
shear terms. In case D3 resistive diffusion is roughly 
three times larger than the sum of the two fluctuating 
terms, while in case D3b the roles are reversed and re- 
solved turbulent dissipation does most of the destruction 
of £tme while the unresolved turbulent dissipation rep- 
resented by our explicit resistivity is relegated to a less 
prominent role. Supporting this transition from unre- 
solved to resolved dissipative processes, in case D3 the 
sum of the fluctuating terms does not show noticeable 
changes in behavior when the magnetic energy is high 
versus when it is low. Instead the response is seen pri- 
marily in the resistive dissipation term. In case D3b, 
however, the fluctuating terms show strong variations in 
response to changes in the magnetic energy. 

This transition from wreath-building dynamos that 
rely on our SGS diffusion to wreath-building dynamos 
that are sufficiently turbulent to be dominated by re- 
solved turbulent dissipation answers one important ques- 
tion relative to the extension of this dynamo mechanism 
to even more turbulent states. It had long been pos- 
tulated that global-scale magnetic structures could not 
exist in the convection zone as they would be quickly 
destroyed by the intense turbulence. While it is clearly 
possible that our wreaths may not be able to survive if 
we were able to simulate far more turbulent conditions, 
case D3b marks an important milestone along the path 
towards the possibility of magnetic wreaths coexisting 
with highly turbulent convection. 

Returning to Figure [TU1 let us now look at the time- 
variation in the production terms. Both cases show vari- 
ability of £tme, but for case D3b we have chosen to show 
a time period that includes states before, during, and af- 
ter a reversal in global magnetic polarity. For both cases, 
careful examination shows that changes in £tme are initi- 
ated primarily by changes in Gms, not by changes to the 
terms dissipating energy. The terms representing both 
resolved and unresolved diffusion respond to changes in 
£tme rather than drive them. In case D3 this is sup- 
ported by the cross-correlation of Gms and Grd peak- 
ing at a 39 day lag, while there is no significant cross- 
correlation between Gms an d either Gfs or Gfa for any 
shift in time. In case D3b both the cross-correlation of 
Gms with Gfs, and Gms with Gfa both peak at a lag 
of 11 days. Resistive diffusion responds faster in case 
D3b with a peak in cross-correlation for a lag of only 
5.6 days. This demonstrates that the variability in the 
toroidal fields is driven by changes in the generation of 
field by the f2-effect. 

If we more closely examine the structure of Gms from 
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Figure 11. Time-evolution of {Aa) between ±45° latitude for case D3 over about 500 days. Times t\ and for case D3 correspond to 
times indicated in Figure [3] Shown are (A^) at (a) the beginning and (b) end of the time interval, (c) the net change between those times 
A(A^), the changes in (A^) due to (d) the fluctuating EMF (Ayl^) pB , (e) the mean EMF (AA^) ME , and (f) resistive diffusion (AA^) RD . 
Of particular importance is the region of positive production in (AA^) FB , which if left unimpeded by diffusion would lead to a reversal in 
global magnetic polarity. The color table has been chosen with a sharp transition from light blue to yellow around zero, thus low-amplitude 
signals, such as seen in (c) and (e), are highlighted. 




Figure 12. Same Figure [TT] but for case D3b. Times t2 and te for case D3b correspond to times indicated in Figures l3l and ITol The 
turbulent EMF induces field of the opposite sense to that which was present at t2 and is opposed by the resistive diffusion. Note that 
(AA^) FE and (A7i0) RD for both cases are topologically similar, but that (AA0) RD is smaller in case D3b, rendering it unable to prevent 
the reversal of (A^) by the fluctuating EMF which begins with the positive region near the equator. 



Equation ([7]). we can expand it to 

d{v+) , {B^)(B e }d(v 4> ) 



Gms 



B^){B r ) 

v v " r r 

r r tan 8 



86 

dV. (13) 



The third and fourth terms are geometric terms from the 
spherical coordinate system which are generally small. In 
order to produce a change in Gms, the dynamo can ei- 
ther change the axisymmetric poloidal field or modify the 
differential rotation of the domain. We have examined 
both the amplitude and geometry of the mean shear due 
to differential rotation and find only very small changes 
in any of the cases presented here. Additionally, reversals 
in the polarity of the wreaths such as those seen in cases 
D3a and D3b require a change in sign for the generation 
term (obtained by dividing by (B^)) and there is never 
a change of sign in the shear profile of the differential 
rotation observed in any of these cases. Thus we are left 
with the conclusion that reversals in the polarity of the 
axisymmetric toroidal fields must be initiated by changes 



in the axisymmetric poloidal fields. 

6.2. Collapse of Resistive Balance Leading to Reversals 

The key to understanding the reversals seen in cases 
D3a and D3b lies in the generation of poloidal field. 
When the poloidal field reverses sign the f2-effect can 
then build wreaths of the opposite polarity and reverse 
the sign of the axisymmetric toroidal field. It is difficult 
to identify a simple model for the generation of poloidal 
field in these cases, particularly in case D3b. We can, 
however, identify the change in the generation mecha- 
nism that occurred betw een cases D3 and D3 b. 

Following the work of iBrown et~aT1 (|2010l |2011[ ). we 
choose to examine the evolution of the <fr component of 
the mean vector magnetic potential (A). This is conve- 
nient as (Af) completely specifies the components of the 
axisymmetric poloidal magnetic field by 



Vx ((A^M =(B r )f+(B e )d. 



(14) 



The temporal variations in the magnetic wreaths arc 
driven by changes in the shear of mean poloidal mag- 
netic fields by mean differential rotation and that only 
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the axisymmctric poloidal fields can change sign, hence 
changes in the polarity of the wreaths can be traced back 
to the evolution of {A$). Further, the key region of the 
domain in which we should monitor (A^) is near the 
equator where the gradients in differential rotation are 
largest and where the wreaths are primarily generated. 
The evolution of (A^) is governed by 

^ = (>> x (B))^ + ((v' x B% - „<J,). (15) 

We have ignored a gauge term in Equation (|15|) which is 
permissible for any longitudinally-periodic gauge. We 
take a time-integral of this equation to look at the 
changes in (A^) over about 500 days in cases D3 and 
D3b, and define the time-integral of each term as 

(AA^) ME = f ((v) x (B)) dt (16) 
(AA 4> ) FE = f\{v' xB'))^dt (17) 

Jtx 

(aa ) rd = - f 2 niJ^dt. (is) 

Thus the change in {A^) can be written as 

A(A^) = (AA^) FE + (AA ) ME + (AA^) RD . (19) 

Figures [TT1 and [T2l show the evolution of (A^) in cases 
D3 and D3b, respectively, as well as the time-integrated 
production of terms shown above and the net change over 
the time interval. For case D3b we chose a time period 
spanning a reversal in global magnetic polarity. In both 
cases (AAj,) ME is small and the evolution is primarily 
governed by the balance between fluctuating EMF and 
resistive diffusion. The primary difference between cases 
D3 and D3b is the collapse of the resistive balance. Both 
cases show similar patterns in (AA^,) FF , namely that the 
fluctuating EMF in both cases is seeking to create a re- 
gion of opposite polarity poloidal field near the equa- 
tor while reinforcing the current sense of poloidal field 
at mid-latitudes. Thus in both cases D3 and D3b the 
turbulent correlations between the existing field and the 
convective turbulence tends to build poloidal field near 
the equator of the opposite sense than the field that built 
the current wreaths through the f2-effect. The difference 
between cases D3 and D3b is that in case D3 the diffu- 
sion term is sufficiently large to prevent the reversal by 
diffusing away the opposite polarity poloidal field at the 
equator before it can accumulate sufficiently to cause a 
reversal. 

What causes the fluctuating EMF to display this 
propensity towards reversing the polarity of (A^) near 
the equator? It would seem that there should be some 
link back to the strong toroidal wreaths, however when 
we expand the fluctuating EMF, we find that 

({v'xB% = (v' r B' d -v' e B' r }. (20) 

Clearly, neither the axisymmetric nor fluctuating com- 
ponents of B,j, come into play here, indicating that 
to complete a reversal we need to connect the large- 
scale toroidal fields to correlations between small-scale 
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Figure 13. (a-c) Values for the three components of a tensor 
relevant to the generation of (e^) as a function of radius and lat- 
itude. Values are computed using a singular value decomposition 
over approximately 3000 days of simulation time with the assump- 
tion that these components of a are spatially local and do not vary 
in time. The a^ r component is very small, whereas the a^g and 
c*0<a components show significant spatial variability and compara- 
ble amplitude, (d-f) Values for components of a t j,j(Bj), showing 
the effect of each component on (e'^}. Magnetic fields have been 
averaged over the same interval as in Figure 1121 (about 480 days). 
Here the contribution of (^^^(B^) is dominant, with a smaller but 
still significant contribution by a^e(Bg). 



poloidal fields and poloidal flows. As shown in Fig- 
ure [TJi), the wreaths are not purely toroidal structures, 
thus the small-scale fields needed in Equation [20] may be 
supplied by the wreaths themselves. However, we have 
not been able to definitively link the poloidal components 
of the wreaths to the reversal process. While the subtle 
nature of this process remains difficult to pin down, we 
do have some hints at its origin. 

6.3. Exploring An a-Like Effect 

The final step in the reversal process is what is often 
described in the pa rlance of mean-field dynamo theory 
as the a-effect fsee ICharbonneaul 120101 ). Generally, the 
a-effect is the source of the axisymmetric component of 
the turbulent EMF, defined as 

(e 7 ) = (v' x B'). (21) 

Specifically, we are interested in the zonal component 
which generates the mean poloidal field and its connec- 
tion to the axisymmetric toroidal field, which might be 
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Figure 14. Magnitude of cross-correlation in time of {e^} and 
{B r } (green), {Bg} (red), and {B^} (blue) for case D3b. Cross- 
correlation is computed as a function of the temporal offset ta,/ 
with negative offsets indicating magnetic fields precede the toroidal 
EMF. Also shown are the 2<r confidence lev els (dashed), compute d 
using a Markov chain Monte Carlo method (Wall & Jenkins 2003). 
The only statistically significant peaks are those relating {e^,} and 



expressed as 

(e' ) = a 0r (B r > + a (j}0 {Bg) + a^B^). (22) 

In its simplest formulation, the components of Q!jj in 
Equation (|22j) are constants, however more complex for- 
mulations exist. 

For case D3b, we have computed the value of the 
three components of the a tensor in Equation (|22[) us- 
ing a singular value decomposition following the work of 
I Racine et al.l (|2011l) . We compute values for a^ at each 
radial and latitudinal location, assuming that ctij is con- 
stant in time. The results of Figure ITBI demonstrate that 
the a ( p c j > (B^) is the most important term in the generation 
of the fluctuating toroidal EMF. Thus, it is particularly 
intriguing to focus on the connection 



(e4) = a^B^). 



(23) 



IBrown et al.l (|2010D showed that for one formulation of 
an a-effect in case D3, was spatially nonlocal, which 
would not be picked up in our fitting procedure. 

The exact mechanism for connecting mean fields and 
the fluctuating EMF is subtle, but we find that in case 
D3b an a-like effect emerges, which is nonlocal in time, 
acting on the same time scale as convective overturning. 
If we consider correlations between the volume-averaged 
magnetic field components and similarity the fluctuating 
toroidal EMF, we find evidence that the a-like effect in 
case D3b is not instantaneous but rather acts on a time 
scale (47 days) which is commensurate with the convec- 
tive overturning time. The volume averages, denoted 
by curly braces, are computed separately for each hemi- 
sphere over all depths and longitudes, and between the 
equator and ±30° in latitude. Combining the data for 
both hemispheres, the cross-correlation is computed and 
shown in Figure [U as a function of the temporal interval 
A T by which {e'^} is offset relative to {B r }, {Bg}, and 
{Be/,} in turn. The peaks in the cross-correlation which 
exceed 2a in significance occur when {e'^} leads {-B^} by 
312 days and when {B^} leads {e^} by 47 days. 

Analysis of the autocorrelation of both {B,p} and{e^} 
indicates that the two peaks are not due to periodicities 
in either of the two time series individually. Further, the 
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Figure 15. Companion plots of the time evolution in case D3b 
of (a) the spherical harmonic coefficients for antisymmetric modes 
with 1 < I < 29 and m = for B^, and (b) {B^) in physical space 
as a function of latitude, both at mid-co nvec tion zon e. Dashed 
lines show times referenced in Figures [3] [4] 1111 and 1121 A factor of 
(— 1)^ — 1 ' / ' 2 is applied to the spherical harmonic coefficients to re- 
move the effect of the wreaths confinement to low latitudes. There 
is a clear progressive spectral transfer of magnetic energy from high 
£ modes to low I modes as each cycle progresses. Reversals begin 
at moderate scales (high I) and then progress to large scales (low 
I). 



widths of these peaks largely originates from the coher- 
ence time for (B^) of about 100 days. The first peak at 
312 days represents the time scale for the Sl-cffect and 
agrees well with the estimate from mean-field theory 
given by 

m ^ p ^MijB^-y (24) 

where P rot is the rotation period, fio is the frame rota- 
tion rate, Af2 is the differential rotation rate, and (B po \) 
is the strength of the poloidal field. For case D3b, this 
yields a value of 324 days. The second peak in the cross- 
correlation between the two fields occurs when ta = —47 
days. This peak suggests that the correlations which gen- 
erate the turbulent zonal EMF are related in some way 
to the axisymmetric toroidal fields, and that whatever 
mechanism establishes this correlation, it has a time scale 
of about 50 days. This temporally and spatially nonlocal 
a-effect clearly points to convection as a key player, as 
other mechanisms like the meridional circulation are at 
least an order of magnitude slower. 

In addition to a timcscale for an a-like effect which is 
commensurate with the convective overturning time, we 
also find evidence for an upscale transfer of magnetic en- 
ergy related to magnetic reversals. Figure ITS] shows the 
temporal evolution of {B^) at mid-convection zone over 
both latitude and spherical harmonic degree i. Several 
reversals of global magnetic polarity arc evident in Fig- 
ure |T5fb) , including the reversal shown in detail in Fig- 
ure|H FigurefTBTa) shows the coefficients of the spherical 
harmonic expansion of the axisymmetric toroidal field for 
the anti-symmetric (odd £) modes with 1 < I < 29 over 
roughly three magnetic activity cycles. In both physical 
space and spectral space, it is clear that each cycle has 
opposite polarity from the preceding cycle. 

There is a preference for antisymmetric modes with 
odd values of £, as would be expected from the fi-effect 
acting on a poloidal field that is preferentially symmet- 
ric about the equator (even I). The upscale cascade 
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Figure 16. Schematic description of the reversal mechanism for cyclic convective dynamos in four steps, (a) Two toroidal wreaths at low 
latitude which generates a turbulent EMF via a nonlocal "a"-effect, either through nonlinear interactions across the equator or via helical 
convection. The sign of the EMF changes at roughly the location of the wreaths, (b) Correlations in turbulent poloidal velocities and 
fluctuating magnetic field drive an induction of mean poloidal field which is roughly octopolar. (c) Mean poloidal field near the equator 
is sheared by differential rotation to generate mean toroidal field through the f2-effect. In these simulations, the largest component is the 
shearing of radial field lines by radial gradients in the differential rotation, (d) Toroidal wreaths of opposite polarity are generated. 



involving odd modes is expected from both theoretical 
and observational studies families of dynamo m odes (see 
iNishikawa fe Kusandl200l iDeRosa et al.ll2012l ). As a re- 
versal occurs we see power showing up first at moderate 
t and then cascading upscale to smaller t values until it 
peaks at I = 3 or 1 depending on the cycle. The reversal 
starts at 25 < I < 29 and then each successive mode 
reverses. There is considerable overlap between cycles, 
in some cases reversals are seen in the high-^ modes in 
as little as a hundred days after the previous reversal is 
completed at \aw-t. We note that convective power peaks 
at spherical harmonic degrees between about 25 and 40 
in these simulations. This suggests that the reversals 
are caused by turbulent processes interacting with the 
wreaths, yielding an upscale energy transfer which or- 
ganizes the large-scale fields. Combined with our cross- 
correlation analysis, this upscale transfer indicates the 
key role of convection in connecting mean toroidal mag- 
netic fields with the fluctuating toroidal EMF. 

As illustrated schematically in Figure [16j the reversal 
mechanism involves three main processes. First, axisym- 
metric wreaths of toroidal magnetic field (Figure ITST a)) 
lead to correlations in the non-axisymmetric poloidal ve- 
locity and magnetic fields which drive an axisymmetric 
turbulent EMF through an a-like effect. The upscale 
transfer of magnetic energy and the fact that the corre- 
lation between the magnetic energy of the wreaths and 
the turbulent EMF peaks on roughly a convective over- 
turning time would seem to point towards the convective 
motions as a key player in this a-like process. In the 
second step, the turbulent EMF reinforces the dominant 
poloidal field at mid-latitudes but is the opposite sign 
near the equator (Figure [16jb)), creating an octopolar 
configuration, with strong radial field concentrations at 
low latitudes (Figure fTBT c^. As the reversal progresses, 
the region of new poloidal field shown in red in Fig- 
ure llGf c) will expand and eventually replace the old sense 
of field shown in blue. The third step involves axisym- 
metric poloidal magnetic field being sheared by differ- 
ential rotation. Here the differential rotation is largely 
cylindrical, thus radial poloidal field is primarily con- 
verted into toroidal magnetic field through the i7-effcct, 



which results in axisymmetric toroidal fields of the op- 
posite polarity (Figure [T6j(d)). The process then repeats 
with the opposite polarity. 

7. TURBULENCE-REGULATED FLUX EMERGENCE 

Photospheric active regions are thought to arise from 
the buoyant destabilization, rise, and emergence of co- 
herent, subsurface toroidal flux structures. It is often 
argued that these subsurface flux structures originate 
below the convection zone, where the strong shear of 
the tachocline promotes toroidal flux generation and the 
subadiabatic stratification of the overshoot region pro- 
motes flux storage by inhibiting the buoyancy in stability 
(jGallowav fc Weisslll98lt Ivan Ballegooiienlll982[ ). In this 
section we offer an alternative viewpoint that is inspired 
and supported by the numerical models presented here. 
Namely, we argue that buoyant flux structures may be 
produced in the Sun and stars not only in the tachocline 
but also in the lower convection zone through the com- 
bined action of rotational shear and turbulent intermit- 
tency. 

In previous papers we have demonstrated that orga- 
nized systems of toroidal flux can persist within a tur- 
bulent convection zo ne despite the inhibiting i nfluence of 
turbulent dispersal (|Brown et al.ll2010l \20Ui} . Here we 
have demonstrated that this continues to hold as we de- 
crease the diffusion, crossing a threshold beyond which 
resolved motions replace artificial dissipation in the dy- 
namical balances that sustain mean flows and fields. Fur- 
thermore, as the diffusion is decreased, intense, localized 
wreath cores form where the magnetic energy density ex- 
ceeds the surrounding kinetic energy density (Figure [2]) . 
This trend is highlighted most dramatically by case S3, 
where the much lower diffusion promotes coherent wreath 
cores stro ng enough to becom e buoyant, as first demon- 
strated bv lNelson et al.l (|2011f ). 

Here we explore in more general terms the link between 
magnetic wreaths and flux emergence, addressing in par- 
ticular on how it might operate in real stars where the 
dissipation is many orders of magnitude less than in sim- 
ulations. We begin by noting that the Sl-cffcct docs not 
just operate on axisymmetric fields; poloidal fields of all 
longitudinal wavenumbers (m) in the convection zone are 
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converted to toroidal fields (of the same m) and ampli- 
fied by rotational shear, blurring the distinction between 
mean and fluctuating fields. Turbulent intcrmittency in 
the surrounding convection can further amplify shear- 
generated flux structures, promoting the generation of 
fibril magnetic fields and coherent, localized wreath cores 
(Figure GO). 

The low Mach number of stellar convection zones en- 
sures that the gas pressure adjusts rapidly to any imbal- 
ance of mechanical and magnetic stresses. Thus, the for- 
mation of fibril, intermittent flux concentrations (wreath 
cores) will induce a pressure perturbation SP ~ Pt — P m , 
where Pt is the turbulent (kinetic plus magnetic) pres- 
sure of the surrounding medium and P m is the mag- 
netic pressure associated with the coherent flux that de- 
fines the wreath core. We have neglected the turbu- 
lent pressure within the wreath core which may be sup- 
pressed by magnetic tension, providing a positive feed- 
back mechanism that can further promote the formation 
of cohere nt, superequiparition wreath cores and buoy- 
ant loops (iKleeorin et alj|198 9; Rog achevskii fe Kleeorinl 
[2007t IK&pvla et al.1120121) . 

Weak magnetic flux concentrations, P m < P t , are not 
susceptible to buoyancy instabilities because their mag- 
netic pressure is insufficient to balance the surrounding 
turbulent pressure, resulting in SP > 0. It is only the 
strongest wreath cores that develop a pressure deficit 
SP < 0, in particular only those cores in which the mag- 
netic pressure P m exceeds the stabilizing influence of the 
surrounding convective motions. This implies that a nec- 
essary but not sufficient condition for the wreath cores to 
become buoyant is that they must be superequipartition 
relative to the surrounding convection. The surround- 
ing flows may in turn enhance or retard the tendency 
for such structures to rise. As demonstrated in Figure [2] 
this is indeed achieved in our simulations and it becomes 
more pronounced as the artificial diffusion is reduced, 
eventually inducing buoyant rise. 

If these superequipartition wreath cores form adiabat- 
ically, this pressure deficit will be accompanied by a den- 
sity deficit e = Sp/p ~ SP/ (7P), established by diverging 
flows along the axis of the wreath core. Radiative heating 
can further warm and rarify the wreath cores, enhancing 
the the density deficit to e ~ SP/P = (P m -P t )/Foaa 
time scale of 

^ - ?ro {'"^'T) < 25 > 

wher e n r is the radiative diffusivity (|Fan fc Fisherl 
119961) . Inserting values from Model S 

(jChristensen-Dalsgaard et al.1 119961 ) for e ~ 1CT 6 
yields r r < 100 days through most of the solar convec- 
tion zone. This value of e corresponds to an emergence 
time r e ~ (2D/eg), of about 10-15 days, where D is 
the depth of the convection zone, and a magnetic field 
strength of B ~ (SireP) 1 ^ 2 - 20-40 kG over and above 
the equipartition value. 

Convection can also promote the buoyant rise of a 
wreath segment by introducing a finite-amplitude un- 
dular displacement, resulting in a draining of fluid 
from from the ap e x of the loop dJouve" fc Brunl [200l 
INelson et all [20111: iWeber et al.l 120111 ). This could in 
principle operate for any field strength but in practice 
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Figure 17. Three-dimensional volume renderings of isosurfaces of 
magnetic field amplitude in case S3. Blue surfaces have amplitudes 
of 10 kG, green surfaces represent 25 kG, and red surfaces indicate 
40 kG fields. Grid lines indicate latitude and longitude at 0.72 
Rq as they would appear from the vantage point of the viewer. 
Small portions of the cores of these wreaths have been amplified to 
field strengths in excess of 40 kG while the majority of the wreaths 
exhibit fields of about 10 kG or roughly in equipartition with the 
mean kinetic energy density (see Figure [2{. 

weak fields will by shredded and reproc essed by convec- 
tion before they emerge (e.g.. lFanll2009( ). 

The dynamics discussed here are indeed exhibited by 
our most turbulent simulation, case S3. Relative to more 
diffusive simulations, this case generates more regions of 
strong, superequipartition fields, as demonstrated in Fig- 
ure[2j and these regions are located in coherent, intermit- 
tent wreath cores, as illustrated in Figure [17] Figure [18] 
highlights two examples in which such wre ath cores be- 
come buoy ant and rise. As d iscussed in INelson et al.l 
(|201lD and INelson et al.l (|2013l ). the loops rise through 
the convection zone through the combined influence of 
magnetic buoyancy and advection, reaching as high as 
0.94i? before they are dissipated by diffusion. The 
wreath which formed these two loops (and two others 
not shown) is not axisymmetric; rather, it spans about 
180° in longitude, reaching peak field strengths of 45 kG. 
We expect the process to be even more efficient in stars 
where the intermittency is presumably much more ex- 
treme. 

In summary, this paradigm of turbulence-induced flux 
emergence postulates that the combined action of turbu- 
lent intermittency and rotational shear generates a broad 
distribution of toroidal magnetic structures and it is only 
the most extreme events, in the high-_B tail of the pdf, 
that become buoyant. It is analogous to the theory of 
turbulence-regulated star formation, whereby supersonic 
turbulence in interstellar molecular clouds generates a 
spectrum of density fluctuations but only the extreme 
events on the tail of the pdf are dense enough to trig- 
ger the Jeans instability an d condense to form protostars 
(jKrumholz fc M cKcc 2005). It is also closely related to 
the negative magnet ic pressur e instability described by 
Kleeorin et al.l (119891) (see also iRogachevskii fc Kleeorinl 
2007tlKapyla et alJl2012tlKemel et al.ll2012D . although it 
does not necessarily rely on the assumptions that under- 
lie that instability analysis, namely scale separation, the 
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Figure 18. Buoyant magnetic loops evolving from small-scale wreath sections amplified by turbulent intermittency. (a) Field line rendering 
of magnetic wreaths at low latitudes in case S3. Field lines are colored by B^ (negative in blue, positive in red) to highlight the two wreaths 
present, (b) Zoom-in on region indicated in (a) showing field line tracings of the core of the buoyant magnetic loops at the same instant 
colored by magnitude of B (weak fields in purple, intense fields in yellow). Volume rendering shows B^ using the same color scheme as in 
(a), (c) The same region 4 days later, showing the continued rise of the loops through the stratified domain and their expansion. 



invariance of the small-scale turbulent energy, and the 
proportionality between variations in the mean and tur- 
bulent magnetic energy (attributed to kinematic shred- 
ding). 

The radial location of the flux bundles that ultimately 
form active regions depends on the kinetic energy den- 
sity in the convection (FKE) relative to that in the dif- 
ferential rotation (DRKE), as well as the efficiency of 
magnetic pumping. In the simulations presented here, 
DRKE/FKE > 1, suggesting that the generation of the 
wreaths is efficient enough that they can persist in the 
convection zone despite magnetic pumping. If this ratio 
falls much below unity, as might be expected for lower 
rotation rates, the wreaths may get pushed toward the 
base of the convection zone. Likewise, if the simulations 
are over or under-estimating the efficiency of magnetic 
pumping, this will influence the location of flux genera- 
tion and the threshold to trigger the magnetic buoyancy 
instability. However, the basic paradigm should still be 
valid. 

The scenario outlined here may resolve several cur- 
rent observational and theoretical puzzles. In par- 
ticular, the non-axisymmetric nature of turbulence- 
induced flux emergen c e is co nsistent with the results of 
IStenflo fe Kosovich ev (2012) who find that many large 
bipolar active regions on the Sun violate Halcs's polar- 
ity rules, and furthermore, that the anti-Hale regions 
often occur at the same latitude as bipoles that obey 
Hale's rules. The fraction of anti-Hale magnetic re- 
gions increases from about 4% for the largest active 
regions (flux $ > 10 23 Mx) to more than 25% for 
smaller bipoles with $ ~ 10 20 Mx. The result that 
more than 70% of intermediate-sized bipoles (4> ~ 10 20 
Mx) obey Hale's laws suggests the presence of organized 
toroidal flux systems throughout the convection zone 
since all of these regions are unlikely to be anchored 
in the tachocline. Meanwhile, the diminishing of mag- 
netic activity patterns with decreasing flux, including an 
increasing fraction of anti-Hale bipoles as well as an in- 
creased scatter in tilt angles and emergence latitudes, is 
often attri buted to the influence of convection on rising 
flux tubes (jJouve fe Brun|[2009t IWeber et al.1l20TTl 12012: 



IJouve et al.ll20Tl . We propose that this intimate cou- 
pling between flux tubes and convection exists not only 
in their rise, but also in their very formation. Finally, 
the non-axisymmetric nature of turbulence-induced flux 
emergence may also account for th e phenomenon of ac- 
tive longitudes (jNelson et al.ll2013l ). 

The observed tilt angles and emergence latitudes of 
bipolar magnetic regions on the Sun is best reproduced 
by models of rising fl ux tubes wi t h initial field str e ngths 
of 20-100kG (e.g. iFanl [200l |j ouve fe Brunl [20091: 
IWeber et ail 1201 it iPinto et all 120111) . However, gener- 
ating such super-cquipartition fields is not a trivial mat- 
ter and in fact represents a formida ble, unresolved prob- 
lem i n solar dynamo theory (e.g., iRempel fe Schusslerl 
1 2 If ) . Laminar amplification of toroidal fields by ro- 
tational shear, the 51-cffcct, tends to saturate at field 
strengths well below cquipar tition due to the back- 
reaction of the Lorentz force (|Vasil fe B rummell 120091 : 
IGuerrero fe KapyTall201lD . Turbulent intermittency can 
help by tapping the energy in the convection that is ulti- 
mately provided by the solar luminosity. It is clear from 
Figure [2] that the coupled action of turbulence and shear 
can generate superequipartition fields of the required am- 
plitude. 

The paradigm proposed here may also help address 
other diffi culties with t achoc line-based dynamos dis- 
cussed by iBrandenburd (|2005D . For example, toroidal 
flux generation does not rely on the radial shear of the 
tachocline, which is maximum near the poles. Instead, 
the expected location of flux generation is where |Vf2 
is maximum in the convection zone. This corresponds 
to the latitudinal shear at mid-latitudes, precisely where 
active regions fir st emerge at the beginning of a cycle, 
as emphasized bv lSpruit] ( 20101) . Note that the potential 
role flux emergence plays in establishing the solar cycle 
is a separate question that we do not address here. 

8. RICHNESS OF STELLAR DYNAMOS 

In this paper we have explored the complex behavior 
of a class of numerical simulations of convective dynamo 
action in rapidly rotating solar-like stars. More broadly, 
however, we have also touched upon the rich landscape of 
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convective dynamo simulations by discussing both per- 
sistent wreath-building dynamos such as cases D3 and 
D3-pml, and cyclic wreath-building dynamos including 
cases D3a, D3b, D3-pm2, and S3. Although the simula- 
tions considered here are ostensibly rotating three times 
faster than the Sun (30© ), the Sun may actually be in 
a similar Rossby number regime, as noted in §1. Thus 
the results presented here may have some bearing on the 
solar dynamo as well as the dynamos of younger, more 
rapidly-rotating solar analogues. 

We have focused on two open questions that arose out 
of our previous work on wreath-building dynamos. The 
first is "Can magnetic wreaths persist in the highly tur- 
bulent conditions of a stellar convection zone" and the 
second is "What physical mechanisms establish and reg- 
ulate the magnetic cycles we see in our simulations?" We 
have also touched upon a third question that all solar and 
stellar dynamo models must eventually face, and that is 
"How are sunspots and bipolar active regions produced 
from dynamo-generated magnetic fields?" . 

The principal issue with regard to the first question 
is whether magnetic wreaths can persist in stellar con- 
vection zones where the magnetic, viscous, and thermal 
diffusion coefficients are many orders of magnitude lower 
than in our simulations. We have investigated this ques- 
tion by systematically decreasing the diffusion in our sim- 
ulations along two complementary paths, one in which 
only the magnetic diffusion coefficient, rj, was reduced, 
and one in which the magnetic and viscous diffusivitiy, 
rj and v, were reduced together, keeping the magnetic 
Prandtl number constant (at a value of 0.5). In both 
cases magnetic wreaths with quasi-cyclic polarity rever- 
sals were attained, although the constant Pm branch ex- 
hibited more regular spatial and temporal behavior and 
thus became the focus of our analysis (see Figure [3]) . 

Although no simulation can approach the extreme pa- 
rameter regimes of stellar interiors, we have demon- 
strated a shift in the dynamical balances that bodes well 
for the possible persistence of magnetic wreaths at much 
higher Reynolds and magnetic Reynolds numbers. In 
short, our simulations suggest that the answer to the first 
question may be "Yes, magnetic wreaths may indeed oc- 
cur in actual stars" . We have investigated in particular 
the balance of angular momentum transport which main- 
tains the differential rotation in our simulations and the 
balance of processes responsible for creating and destroy- 
ing the magnetic energy of the wreaths. In both cases, 
as we move from case D3 to case D3b we find that re- 
solved turbulent dissipation has taken the place of SGS 
dissipation (see Figures 181 and fTUj) . This is an important 
milestone towards demonstrating that wreaths can exist 
in highly turbulent settings and that they are not reliant 
on the explicit diffusion in previous simulations. 

We have found that magnetic wreaths persist in our 
higher-resolution, lower-dissipation, more turbulent sim- 
ulations, yet their nature is altered in a fundamental and 
significant way. Most notably, they are no longer axisym- 
metric. In our more turbulent simulations such as case 
D3b, the nearly axisymmetric wreaths of case D3 are re- 
placed by coherent wreath segments, typically spanning 
between 45° and 270° in longitude. This is associated 
with a shift in the magnetic power spectrum from lon- 
gitudinal wavenumber to = to moderate to values. It 
also has important implications for flux emergence, as 



discussed with regard to question 3 below. 

The first clues as to the origin of magnetic cycles in our 
simulations (question 2) were uncovered by iBrown et al.l 
(|2010L 120111 ). showing that one can move from a per- 
sistent wreath-building dynamo state to a cyclic one by 
increasing the rotation rate. Here we have shown that a 
similar transition from persistent to cyclic wreaths can be 
achieved by decreasing the effective magnetic diffusion, 
and thereby increasing the magnetic Reynolds number at 
a fixed rotation rate. As mentioned above, the constant- 
Pm branch of solutions exhibited more regular cyclic be- 
havior despite the higher degree of turbulence. 

We have not obtained a definitive exposition of the 
physical mechanisms that give rise to and regulate mag- 
netic reversals. However, we have traced their opera- 
tion to the zonal component of the turbulent electromo- 
tive force (EMF) near the equator. In case D3 diffusion 
prevented reversals in the polarity of the axisymmetric 
poloidal field by locally offsetting the creation of mean 
poloidal field by turbulent fluctuations. In the lower- 
dissipation case D3b, this balance breaks down, leaving 
a residual turbulent EMF near the equator that creates 
poloidal field with a polarity that is opposite to that 
of the pre-existing field, as shown in Figure [12l Once 
magnetic reversals are thus initiated, the overall reversal 
process follows the schematic description found in Fig- 
ure [16] 

Our simulations cannot directly address the third ques- 
tion regarding how solar and stellar dynamos produce 
sunspots and bipolar active regions. The detailed dy- 
namics of flux emergence are too intricate to reliably cap- 
ture in any current global dynamo simulation. However, 
the change in the nature of the wreaths as the dynami- 
cal balances shift suggests that they may play an impor- 
tant role in generating buoyant magnetic loops in actual 
stars. As discussed in §7, these simulations suggest that 
strong, coherent magnetic structures of moderate angu- 
lar extent can be created in the cores of the magnetic 
wreaths. If this trend were to continue to the extremely 
low diffusion regimes of actual stellar convection zones, 
one would expect these flux bundles to become buoyantly 
unstable and rise. Indeed, this expectation is confirmed 
by our simulation Case S3 that employs a less diffusive 
SGS model and that exhibits the self-consistent genera- 
tion of buoyant toroidal flux tubes in a global convective 
dynamo simulation. This picture of flux emergence as a 
fundamentally turbulent process contrasts strongly with 
more idealized scenarios where the principal role of con- 
vection is simply to produce a differential rotation. One 
might expect this revised paradigm to have observable 
consequences in such active region characteristics as dis- 
tribution, tilt angle, and helicity. Furthermore, it may 
call into question our traditional reliance on sunspots as a 
straightforward proxy for the axisymmetric toroidal field 
at or below the base of the convection zone. 

The rich behavior of these systems provides important 
insight into the dynamo models for the Sun and solar- 
type stars. The trend towards non-axisymmetric fields 
with enhanced turbulence, while still maintaining global- 
scale organization, pushes at the boundaries of our un- 
derstanding of dynamo theory in solar-like settings. That 
these mechanisms are accessible with current computa- 
tional resources clearly invites further intensive study of 
these topics. 
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APPENDIX 

A. THE DYNAMIC SMAGORINSKY SGS MODEL IN ASH 

We implement the dynamic Smagorinsky sub grid-scale 
model in the ASH c ode following the prescription of 
iGermano et al.l (119911). inspir ed on the original formu- 
lation by iSmagorinskvl (|1963l ). The Smagorinsky eddy 
viscosity v$ is defined as 



,1/2 



(Al) 



where ey is the resolved strain-rate tensor, A is the 
grid spacing, and Cs is the constant of proportional- 
ity. The dynamic Smagorinsky model (used here in case 
S3) makes an assumption of self-similar behavior in a re- 
solved inertial range of a turbulent cascade to determine 
the value of Cs at each point in the domain. This is 
done by choosing a test-scale which is larger than the 
grid-scale by a factor /3. Traditionally, and in this work, 
[3 = 2. Only high resolution ASH simulations are able to 
assure that /3A is suitably within an inertial range. 

If we define the full velocity field, including scales that 
are not resolved in the simulation, as m, then the hat op- 
erator denotes grid-scale filtering keeping only the scales 
that are resolved in the simulation and the tilde opera- 
tor expresses filtering at the test scale keeping only scales 
that are resolved by the simulation and larger that /3A. 
We define the stress tensor at the grid-scale A as 

Sjj = UiUj — Uiiij (A2) 

and the stress tensor at the test filter scale /3A as 

Sij = UiiTj — Uiilj. (A3) 
Finally, the resolved stresses can be written as 

Sij = UjUj — Uiilj. (A4) 

Note that Sij can be computed directly from the resolved 

flows while Sij and Sij require a SGS model, such as the 
Smagorinsky model. 

By assuming scale-invariance in the inertial range of 
the turbulent spectra, one obtains 



Note that this is simply stating that the resolved stresses 
are the difference between the stresses filtered at the test 
scale and the stresses filtered at the grid scale. Apply- 
ing the Smagorinsky model and contracting with to 
obtain a scalar equation results in 

Sij&ij = ~2C S (/? 2 A 2 
Solving for Cs gives 
C s 



A 



2| e, 



Sij&ij 



2A 2 e l? /3 2 



(A6) 
(A7) 



To assure computational stability we require Cs to be 
positive and apply a spectral filter on the denominator 
of Equation (|A7[) which removes scales smaller than /3A. 
Additionally a Gaussian smoothing with a width equal 
to the largest horizontal grid spacing is applied to Cs 
to prevent grid-scale ringing in the Cs field. With the 
SGS viscosity thus determined, we apply constant SGS 
Prandtl and magnetic Prandtl numbers in order to deter- 
mine the SGS thermal diffusivity and magnetic resistivity 
coefficients at each point in the domain. 

The dynamic procedure involving a turbulent cas- 
cade has known problems adjacent to impenetrable walls 
where viscous boundary layers form that are rather lam- 
inar, such as the upper and lower boundaries in ASH 
(|Meneveau fc K atz 2000(1 . To compensate for a dimin- 
ished cascade, we introduce a smoothly varying enhanced 
eddy viscosity immediately adjacent to there boundaries, 
occupying only 0.3% of the domain in radial extent. 

In the dynamic Smagorinsky model the nonlinear na- 
ture of the diffusion term requires an explicit time step- 
ping scheme which imposes an upper limit on the size 
of our time-step. In order to control the time-step, an 
artificial ceiling is placed on the dynamic Smagorinsky 
viscosity. In case S3, on average, the ceiling is applied 
to 800 out of 76 million grid points at each time step. 
As the time step is required to be less than A 2 lin /f for 
numerical stability, the functional form of the ceiling is 
given by 



t A z ■ 



(A8) 



where A m j n is the smallest local grid-spacing in any di- 
rection and i max is the desired size time step. In case S3 



B. GENERATION OF DIFFERENTIAL ROTATION 
KINETIC ENERGY 

As shown in Equation ([2]). the time evolution of an- 
gular momentum in our domain can be written in con- 
servative form as the divergence of a flux vector JF. The 
radial component is given by 



T r = p\ 



VT ^r (~) + V ' 4>V ' r + + ^ r ^ ^ 



'~a — -BLBr ~ ~. — -B<l>B r 



(Bl) 



Sh 



(A5) 



where the terms are from left to right due to viscous dif- 
fusion, fluctuating Reynolds stress, mean Reynolds stress 
from the meridional circulation, the Coriolis force with 
Oo representing the frame rotation rate, the Maxwell 
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stress, and mean magnetic torques. The latitudinal com- 
ponent is given by 



Fg = pX 



z^sin^ d 
r~dl 



\smOJ 



+ v '<h v 'e + fa v<f, + vgfl X 



47T/3 



-B'B' 



47T/9 



B^Bg 



(B2) 



where the terms have the same ordering and identities 
as in the radial component. We ignore the flux due to 
the Coriolis force because while it can be large locally, 
it cannot do any net work on the system when averaged 
over the full domain. We can also write the fluxes in 
cylindrical coordinates in terms of the cylindrical radius 
A and the distance from the equatorial plane z. The flux 
in cylindrical radius is given by 



F\X = F r sin Of + Fg cos ( 
while the flux in z is given by 



(B3) 



F z z = F r cos Of ~ F e sin 90, (B4) 

If we multiply equation ([2]) by the longitude-averaged 
rotation profile Cl, we are left with an equation for the 
time evolution of the kinetic energy density in the mean 
differential rotation profile (Edr), 



dt 



O V- F 



(B5) 



We take a volume integral over the entire domain in order 
to calculate the total rate of change in the kinetic energy 
of differential rotation and rewrite the right-hand side as 



d(E BR ) 
dt 



dV 



F-Vn-V ■ [CIF 



dV. (B6) 



The second term in the integral can be rewritten using 
the divergence theorem as a surface integral, leaving us 
with 



dt 



dV= / F-VildV- / {lF r dS. (B7) 
Jv Js 

Our choice of impenetrable and stress-free boundaries 
causes all of the hydrodynamic terms terms in the surface 
integral to vanish on both the inner and outer bound- 
aries. Likewise our choice of a perfect conductor bound- 
ary condition on the lower surface causes both the fluc- 
tuating and mean magnetic torques to vanish there. The 
choice of a potential field boundary condition on the up- 
per surface forces the mean magnetic torques to be ex- 
actly zero, however it does in principle allow the Maxwell 
stress to be non-zero. This reduces the surface integral 
to 



flF r dS = 



^B^Rlsm'OdO, 



(B8) 



We have calculated this term to be about five orders 
of magnit ude smaller than the volume integral term in 
Equation (|B7P when averaged over long periods in cases 
D3, D3a, and D3b. We chose to ignore this surface term 
in our analysis of time-averaged quantities. 

The generation and dissipation of differential rotation 
kinetic energy can be written as the sum of five terms, 



as was done in Equation Those terms, which repre- 
sent viscous diffusion, Reynolds stress, meridional circu- 
lations, Maxwell stress, and mean magnetic torques, are 
given in turn by 
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The time-averaged values of these terms are reported in 
Tableland Figure M 
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